Exploring trends in B/Bmsy.
library(here)
## here() starts at /Users/frazier/github/ohiprep_v2021
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.5 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.0.2 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(broom)
key <- read_csv(here("globalprep/fis/v2021/output/fis_bbmsy_gf.csv"))
## Rows: 609501 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): stock_id, bmsy_data_source, RAM_gapfilled
## dbl (5): rgn_id, taxon_key, year, bbmsy, mean_catch
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
scores <- read_csv(here("globalprep/fis/v2021/output/fis_bbmsy.csv"))
## Rows: 146602 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): stock_id
## dbl (3): rgn_id, year, bbmsy
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
catch <- read_csv(here("globalprep/fis/v2021/output/stock_catch.csv"))
## Rows: 90333 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): TaxonName, CommonName, Resilience, stock_id
## dbl (4): year, fao_rgn, TaxonKey, tons
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
bbmsy <- left_join(scores, key) %>%
filter(bmsy_data_source=="RAM") %>%
filter(is.na(RAM_gapfilled))
## Joining, by = c("rgn_id", "stock_id", "year", "bbmsy")
head(bbmsy)
## # A tibble: 6 × 8
## rgn_id stock_id year bbmsy taxon_key bmsy_data_source RAM_gapfilled
## <dbl> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 1 Istiompax_indica-… 2001 1.49 600217 RAM <NA>
## 2 1 Istiompax_indica-… 2002 1.47 600217 RAM <NA>
## 3 1 Istiompax_indica-… 2003 1.47 600217 RAM <NA>
## 4 1 Istiompax_indica-… 2004 1.44 600217 RAM <NA>
## 5 1 Istiompax_indica-… 2005 1.34 600217 RAM <NA>
## 6 1 Istiompax_indica-… 2006 1.35 600217 RAM <NA>
## # … with 1 more variable: mean_catch <dbl>
bbmsy <- bbmsy %>%
select(stock_id, year, bbmsy) %>%
unique()
bbmsy_trends <- bbmsy %>%
group_by(stock_id) %>%
do(mdl = tidy(lm(bbmsy ~ year, data = .))) %>%
unnest(mdl) #fitting a model per each goal and rgn_id
## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
# data.frame(glance(data_lm, mdl))
bbmsy_trends2 <- bbmsy_trends %>%
filter(term == "year") %>%
select(stock_id, term, estimate, p.value) %>%
data.frame() %>%
filter(estimate <= 2) %>%
filter(estimate >= -2)
hist(bbmsy_trends2$estimate)
stocks <- unique(bbmsy$stock_id)
for ( stock in stocks){ # stock <- stocks[1]
f <- ggplot(filter(bbmsy, stock_id==stock), aes(x=year, y = bbmsy)) +
geom_point() +
geom_line() +
xlim(c(2000, 2020)) +
labs(title=stock) +
theme_minimal()
print(f)
}
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?